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Abstract 

We review the methods developed for combining the parton shower approximation 
04 ■ to QCD with fixed-order perturbation theory, in such a way as to achieve next-to- 

leading-order (NLO) accuracy for inclusive observables. This has made it possible to 
generate fully-simulated hadronic final states with the precision and stability of NLO 
calculations. We explain the underlying theory of the existing methods, MCQNLO and 
^ I POWHEG, together with their similarities, differences, achievements and limitations. 

D . For illustration we mainly compare results on Higgs boson production at the LHC, with 

particular emphasis on the residual uncertainties arising from the different treatment 
, of effects beyond NLO. We also briefiy summarize the difference between these NLO 

-I- parton shower methods and matrix-element + parton shower matching, and current 
efforts to combine the two approaches. 
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1 Introduction 



In the past, parton shower generators and next-to-leading-order (NLO) calculations were 
seen as complementary approaches to computing hadronic interactions. The former had 
the more practical purpose of assisting experimental physicists in planning and carrying 
out experimental analysis. The latter were instead aimed at performing precision tests of 
perturbative QCD. The simplicity of the framework of fixed-order calculations was in fact 
required in order to make unbiased comparisons of theoretical calculations and data. After 
the many tests of QCD carried out at lepton and hadron colliders, convincing evidence was 
established that perturbative QCD at the NLO level works well, and improves the agreement 
of theoretical prediction with data. The interest in precise NLO calculations has then shifted 
in the direction of predicting cross sections and backgrounds for collider processes. In the 
meantime, a theoretical effort has begun aimed at improving shower generators with the 
use of more precise matrix elements. This effort has led to two new developments that have 
had a significant impact on collider phenomenology: matrix-element and shower matching 
(ME+PS), and NLO calculations interfaced with showers (NLO-I-PS). The former originated 
from the so-called CKKW paper [T| , and several implementations and variants of the original 
method have subsequently appeared in the literature (see [2] for a review). In the present 
review, we focus upon the NLO+PS development, initiated by the MC@NLO paper [3] and 
followed later by the POWHEG proposal . The aim of these methods was to improve the 
event generation of a basic process in such a way that NLO accuracy is reached for inclusive 
observables, while maintaining the leading logarithmic accuracy of the shower approach. 

This review is organized as follows. In Section [2] we briefly recall the structure of an 
NLO calculation, with emphasis on those aspects that are needed for the implementation of 
NLO-I-PS generators, and in Section[2]we review the basics of parton showers. In Section^] 
we contrast the description of the Sudakov region in fixed-order calculation and in shower 
algorithms, and specify the requirements for a NLO-fPS generator. In these sections, for 
concreteness, we will always make reference to the simple example of Higgs boson production 
in giuon fusion in order to clarify the basic concepts. 

Sections [5] and [6] illustrate the basics of the MC@NLO and POWHEG methods. The 
following three sections discuss few issues that are common to both methods. In Section [7] 
we discuss the use of parton density functions in the NLO-fPS framework, and in Section [SJ 
we discuss a commonly used method for the inclusion of spin correlations. 

The illustration of the methods given in Sections [5] and El are limited for simplicity to 
cases where there is only one singular region to be considered. There are in essence no 
conceptual difficulties in dealing with more complex cases, and in Section [5] we give some 
indications of how this is done. 

Section (TU] illustrates the role of truncated showers, which were introduced in ref. [3] as 
a requirement to preserve colour coherence in the POWHEG approach. In this framework, 
the relation between POWHEG and MC@NLO is also better clarified. This common view 
of both methods allows one to better understand the origin of differences between them, as 
discussed in Section [TT] In Section [T^l the dominant uncertainties in NLO-I-PS generators 
are discussed, using again as an example the process of Higgs production via gluon fusion. 

In Section [T3] we compare the NLO-I-PS and ME-I-PS approaches, with the aim of clari- 
fying when either approach should be preferred. 

Finally, in Section [14] we briefly summarize the future direction of improvement in the 
development of NLO-I-PS generators. 

2 Next-to-leading-order calculations in QCD 

Next-to-leading-order (NLO) calculations in QCD are used to compute infrared and coUinear 
safe quantities at the one loop level. They considerably reduce the uncertainties of theo- 
retical predictions, and experience from e'^e" , ep and hadron colliders has shown that they 
lead to remarkable agreement of theory with data. Since QCD radiation has coUinear and 
infrared divergences, infrared and collincar insensitivity is an unavoidable requirement for 
an observable to be perturbatively calculable. Thus, fixed-order NLO calculations cannot 
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be used for fully exclusive observables, and are not straightforwardly interfaced to parton 
shower programs. 

A detail explanation of the general structure of an NLO calculation, with particular 
attention to its use in the context of shower matching, is given in ref. The reader will 
also find there references to the most popular subtraction methods for the implementation of 
QCD NLO corrections. Here we will give an elementary illustration of the typical structure 
of an NLO calculation, using as an example the gluon fusion production of an on-shell Higgs 
boson in hadronic collisions. In this example, the Born phase space is characterized by a 
single variable that can be taken to be equal to the Higgs rapidity y. We write the Born cross 
section as B d<I>B , where B is the differential Born cross section, and $b is the Born phase 
space. In our Higgs example, B = daH/dy, and d$B = dy. The calculation of the NLO 
cross section for the production process requires the inclusion of the virtual corrections V 
and of a real emission process, with a phase space including one extra parton with respect to 
the basic process. In our example, the real production process is the production of the Higgs 
plus one extra gluon. We neglect at this point qg, gq and qq processes that arise at NLO, 
in order to keep our illustration as simple as possible. The real cross section kinematics can 
be characterized by the Higgs rapidity, as before, and by three variables associated with the 
emission of an extra parton. We can choose for these variables the cosine of the angle of the 
emitted parton momentum with respect to the beam direction, its energy and its azimuth, 
all measured in the Higgs+parton rest frame (i.e. in the partonic centre-of-mass system). 
We write the real cross section as Rd^R, where d^^ = dy dcos6 dE dcfi. We will also write 
d$fl — d$Bd$rad, with d^iad = d COS 6* d£' d0. We call $rad the radiation phase space. The 
differential NLO cross section is written schematically as 

da = (b($b) + 1>($b)) d$B + i?($R)d$R , (1) 

where we have also assumed that all phase space Jacobians needed to reproduce the Lorentz 
invariant phase space are absorbed into the definition of B, V and R. The virtual cross 
section V (which we assume to be renormalized) has soft and coUinear divergences, and 
must thus be regulated with an infrared cutoff, or by using dimensional regularization. The 
real cross section also exhibits infrared and ultraviolet divergences after integration over the 
radiation phase space, such that the full cross section is finite after integration. We have 
not included in Equation [T] a collinear counterterm, which is needed to remove initial-state 
collinear singularities, since our aim here is to illustrate the structure of the NLO corrections, 
rather than give a detail explanation of how they are structured. 

Notice that we assume that the real phase space $r can be given as a function of the 
Born phase space and a radiation phase space, $r = $r($b, 'I'rad)- This is straightforwardly 
achieved in our example, but it can be done in general. We will say that ^b is the underlying 
Born configuration of the real phase space. We also require that the mapping $r($b5 ^'rad) 
has the following properties: in the limit in which the radiated parton is soft, <I>r coincides 
with $B after removal of the soft parton, and in the limit in which the radiated parton is 
collinear to another massless parton, <I>r coincides with $b after merging of the collinear 
partons. By merging we mean that the two collinear partons are replaced by a single 
parton, with momentum equal to the sum of the two in case of final state radiation, or 
equal to the difference of the two in case of initial state radiation. In the Higgs example, 
the underlying Born configuration is obtained by performing a longitudinal boost to the 
frame where the Higgs has zero rapidity, and then performing a transverse boost, such that 
the Higgs transverse momentum vanishes. After that, the inverse of the initial longitudinal 
boost is applied. In this way the Higgs rapidity remains equal to that in the underlying 
Born configuration. 

The cancellation of soft and collinear singularities in NLO calculations is usually dealt 
with using the so-called subtraction method. We introduce an auxiliary counterterm C(<I>r ), 
which is required to coincide with the real squared amplitude R in the soft and collinear 
limits. Assume now that we want to compute an infrared safe observable O. Infrared safety 
requires that 

0($R($B,*rad)) ^0($b) (2) 
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in the soft or collinear limit. We can write 



(o) = y"d$B(s($B)+n$B))o($B)+yd$R,i?($R)o($R) 

= j d$B [B($b) + Vi^B)] 0($b) + J d$R [i?($R)0($R) - C($r)0($b)] (3) 

where 

F($b) = V{^b) + J d$rad C($r($b, $rad)) ■ (4) 

The above equations represent schematicahy the subtraction method in QCD. By a suitable 
choice of the counterterm C, the integral of the radiation variables in Equation |4] can be 
performed analytically. The soft and collinear divergent terms arising from this integration 
cancel against those of the virtual term, V, yielding a finite result V in the sum. At the 
same time, in Equation [3] the soft and collinear divergences in R cancel, because in the soft 
or collinear limit 0($r) = 0(<i>B), and C has the same singularity structure as R. 



3 Parton shower approximation 

The collinear and infrared divergences of QCD are associated with enhanced amplitudes in 
collinear and soft regions of phase space. These enhancements are manifest in jet production, 
and in the rapid increase in particle multiplicity in hard scattering processes. They give 
rise to large logarithmic coefficients in observables that involve widely different scales, for 
example the jet mass M at a hard process scale Q ^ M. The parton shower approximation 
aims to take the enhanced contributions into account to all orders, neglecting terms with 
subleading logarithmic or constant coefficients. 

We consider first the leading collinear region in which an extra parton is emitted at a 
small angle by one of the outgoing lines of an n- parton process (final-state emission). Here 
the cross section approximately factorizes in the form 

dcr„+i($„+i) = V{^ri,d) d<Tn{^n) d$rad (5) 

and d$„-|-i — d$„d$rad- The function V depends on the type of emitting and emitted 
partons. In the notation of the previous Section, if the n-parton process is the Born process, 
we have $„ = $b, = ^r, and 7'($rad) B($b) = i?'^^^'($R) is an approximation to 

i?($R) in the near-collinear region. In this region it is convenient to parametrize the phase 
space $rad in terms of a hardness scale q, for example the transverse momentum relative 
to the emitter, the fraction of longitudinal momentum z of the emitter after the emission, 
and the azimuthal angle of splitting, cj). Then 

Wad)d$,.,«^^P(z,0)dz^, (6) 

n q ZiT 

where P{z, </>) is the relevant (DGLAP) splitting function, which in practice is averaged over 
the azimuthal angle and simply written as P{z). The collinear divergence is regulated with 
a cutoff, q > Qq. Emissions with q < Qq are said to be unresolvable. Emissions with small 
momentum fractions z are also unresolvable. Depending on the definition of the scale q, the 
cutoff on z is some function zo{q, Qo)- 

Equations [S] and IHl provide the basis of an iterative scheme for summing coUinear- 
enhanced contributions to all orders. Each parton participating in a hard process can emit 
at scales q up to order Q. The probability of an emission in an interval between q + dq and 
q is given by Equation [S] It follows that the probability of no resolvable emissions between 
scales qi and (72 < 91 is given by 



As (91,92) = exp 



q2 ^ 1 Jzo 



(7) 



This function is known as the Sudakov form factor. 
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We have justified Equation[7]using unitarity, i.e. the conservation of probability. We now 
remark that relying solely upon field theory would lead to the same conclusions. In other 
words, it can be proven that the inclusion of all logarithmically enhanced virtual corrections 
in a shower process amounts to the inclusion of running couplings at each splitting vertex, 
evaluated at a scale of the order of the virtuality of the incoming parton, supplemented by 
the insertion of a Sudakov form factor on each internal line, having as argument its vir- 
tuality. The result satisfies unitarity, which we do expect, since field theory must satisfy 
unitarity. We can say, more specifically, that logarithmic singularities must cancel in in- 
clusive quantities, thanks to the Kinoshita-Lee-Nauenberg (KLN) theorem [3 [7], which in 
turn is a consequence of unitarity. Thus we expect all singularities to cancel order by order 
in perturbation theory, and in fact we see that by expanding the Sudakov form factor in 
Equation[7]at order as we do get a term that cancels exactly the integral of the real diagram 
describing the splitting process in the collinear approximation, given by Equations [S] and [Bl 

In a parton shower event generator, multiple emissions are generated by the Monte Carlo 
method. Given the initial scale Q, the scale of the first emission is found by solving the 
equation As((5, qi) = Ri, where Ri is a pseudorandom number uniform on the interval [0,1]. 
The next emission is at q2 < qi where As((7i,(i'2) = ^2, and so on, until the chosen scale 
falls below the cutoff Qo . The emitted partons can themselves emit in the same way. In this 
way each outgoing parton from the hard process is converted into a shower of partons at the 
cutoff scale Qo- Below this scale the running coupling becomes large and colour confinement 
ensures that the parton showers are converted into jets of hadrons. Perturbation theory is 
not applicable to this process and one has to resort to hadronization models, which fall 
outside the scope of the present article (for a review see [2]). One important feature of the 
hadronization process, observed experimentally and built into the models, should however 
be noted here, namely its locality: the flow of energy-momentum and flavour in the parton 
showers is preserved, up to power-suppressed corrections, in the resulting hadron jets. This 
justifies the comparison of results in perturbative QCD with data on hadronic flnal states. 

Shower Monte Carlo generators also treat initial-state radiation in the collinear approx- 
imation, but the kinematics are slightly different: a parton from an incoming hadron beam 
starts at a low scale with its momentum along the beam direction and recoils with a more 
spacelike four-momentum after each emission, eventually reaching the hard process scale 
q Q with a significant transverse momentum. For technical reasons it is more convenient 
to treat this in the opposite direction, evolving incoming partons from the hard scale down 
to the hadronic scale, where they are matched to the parton distributions of the incoming 
hadrons. 

In order to arrive at a formula for the backward splitting probability in initial-state 
radiation, we begin by writing the forward splitting probability in the form 

7'($,.ad)d$,ad ~ ;^d$,ad , (8) 

where R^^i and R^'~^ are the n+1 and n body cross section in the MC approximation, in the 
singular limit. This formula can be immediately generalized to initial-state radiation, the 
only difference being due to the fact that the ratio of parton density functions in R^^^/R^'^ 
does not cancel as in final-state radiation, leading to the formula 

P(-«)($.ad)d$,ad « ^^p(,)/^|Zfl^dz^, (9) 

TT q j{x,q) 27r 

which is Sjostrand's formula for backward evolution [Sj in initial state radiation. 

Collinear parton emission can be treated at the cross-section level, as in Equation [5l 
because an emission cannot be collinear with more than one hard parton. This is not the 
case for soft gluon emission, which occurs coherently from different partons and so should be 
treated at the amplitude level. Two approximate showering schemes have been developed 
to take account of soft gluon interference. In angular- ordered showering, the scale variable q 
is defined so that successive parton emissions are at decreasing angles. The separate colour 
charges of a pair of partons in the shower with opening angle 9 can be resolved by soft 
gluons emitted at smaller angles, 0' < 9, so this emission is treated incoherently. On the 
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other hand, when 9' > the partons are not resolved and they emit as a single object with 
the colour charge of their parent parton. Therefore this emission is treated as coming from 
the parent, corresponding to angular ordering of successive emissions, as shown in Figure [TJ 




Figure 1: Angular ordering of coherent soft gluon radiation (6*' > 9). 

The other method for treating soft gluon coherence is dipole showering [9]. Instead of 
sequential splitting of one parton into two, pairs of partons are treated as dipole sources for 
gluon emission, which splits one dipole into two. This is valid in the large- TV approximation, 
where N is the number of colours. Some subleading colour effects are included by matching 
the resulting cross sections to the parton splitting functions in coUinear regions. 



4 The Sudakov region 

The most visible difference between an NLO result and the output of a parton shower is the 
structure of the so-called Sudakov region, i.e. the region where radiation with low energy 
and transverse momentum is important. In the example of Higgs production, this is the 
region where the transverse momentum of the Higgs is much smaller than its mass. The 
©(ttg) result for this quantity diverges as \og{pT / Mh) / Pt in the pr ^ limit. However, 
the integral of the pt distribution from zero to any finite value is finite. In other words, 
one should imagine that a contribution formed by a negative i5 function with an infinite 
coefficient is located at pt equal to zero, and that it cancels the divergent contribution 
arising from the real emission cross section. 

In the small pt region, perturbation theory is not reliable any more, and one needs 
to resum the large logarithms of pt to all orders in perturbation theory. In the leading 
logarithmic approximation, this distribution is dominated by recoil of the Higgs against the 
hardest emission. The shower Monte Carlo algorithm yields for the hardest emission the 
cross section 



, , = -— (5(pt)A(Qo) + A{pt)-— — , A(pt) = exp 
dy dpT dy dy dpT 



Q 



di 



(MC) 



dy 



(10) 



The factor A(pt) arises from the product of the two As(Q,pt) Sudakov form factors, 
associated with the radiation from each initial line, in the notation of Equation [71 Here Q 
is the hard process scale, of the order of A/h hi this case, and cr^'^^-' denotes the shower 
Monte Carlo approximation for the real emission cross section. If we assume that the shower 
algorithm is ordered in px, Equation [TU] is a simple consequence of the shower formula, since 
it corresponds to the generation of the first emission by the algorithm. It can be shown, 
however, that the formula also holds for angular ordered showers, where the hardest emission 
is not necessarily the first [4 . We notice that the shower unitarity relation holds 



r ^(pt)A(Qo) + Zl;^ dpT = A(Qo)+r^d^.T = A(Q) = l(ll) 
^0 ^ Jqo dpT 



so that 







dydpT dy 



This relation is independent of the particular form of cr'-'^'-'^ . 
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A shower Monte Carlo program alone will generate a transverse momentum distribution 
that is accurate only for small transverse momenta, since dcr^^^^ is reliable only in the 
coUinear approximation. For small transverse momenta, however, rather than having the 
singular behaviour of an NLO calculation, it is well behaved, with the Sudakov form factor 
damping the small pt singularity of the tree level result. Many event generators are capable 
of adding a matrix- element correction (MEC), such that for large transverse momentum the 
shower result matches the fixed-order result |10j . This is achieved, in essence, by replacing 

^(NLO) in 

Equation [TUl Assuming, for the moment, that we are dealing with a 
shower algorithm ordered in transverse momentum, the generation of the first emission in 
MEC is given by 



^^(MEC) ^ 



A(0o 



R 

A(pT)-gd$rad 



A(pt) = exp 



rad 



(13) 

The notation used in Equation [13] deserves some explanation. We write in a compact nota- 
tion a fully differential cross section that can have different final states as a single formula. 
The first term in the square bracket represents the production of an event with the Born 
kinematics, and phase space $b- In the Higgs example, it represents a Higgs boson with zero 
transverse momentum. The second term represents the full real process, with production of 
a Higgs and a parton, balanced in transverse momentum. The above formula represents the 
probability that either event is produced. 

The shower unitarity Equation [TT] is then written in the general form 



A(go) 



R 

A(pT)^d$rad = 1 



(14) 



where it is intended that the d$rad integration is limited to the region where pt{^r) > Qo- 
In Figure [5] we give a pictorial representation of the distribution of the transverse mo- 
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Figure 2: Transverse momentum distribution of the Higgs at NLO, in a shower algorithm, 
and in a MEC shower. 

mentuni of the Higgs boson at fixed rapidity at NLO order (i.e. 0{a^)), from the shower 
algorithm, and from a MEC shower algorithm. For the NLO result, one should imagine that 
the NLO curve diverges at small px up to a tiny cutoff, and that a tiny bin with a very 
large, negative value is located at pt = 0. The resummation of coUinear and soft singulari- 
ties performed by the shower algorithm using the exact real emission cross section starts to 
differ from the LO one at pt around 40 GeV, and for smaller pT it tames the divergence of 
the NLO cross section. The shower approximation has the same behaviour for moderate to 
small Pt, but it drops rapidly as pt approaches the maximum scale of radiation allowed by 
the shower algorithm (an exact implementation of Equation [10] would imply that the cross 
section vanishes exactly for pt > Q. Subsequent emissions in the shower process will tend 
to smear the region of ~ Q)- The area under the two shower curves equals the Born 
cross section. 

The main objective of a NLO-I-PS implementation is to improve the shower approxima- 
tion, in such away that it achieves NLO accuracy for inclusive quantities. Thus, referring to 
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Figure[2l we would expect that the Higgs transverse momentum distribution in an NLO+PS 
approach should have the smooth shape of the Shower and MEC approaches at small p^, 
should match (up to even higher order terms) the NLO calculation at large p^, and should 
have the total area of the curve equal to d<j^^'~' / dy . In other words 

, , H ^ K(y,pr^)——^ . (15) 
dy dpj, dy dpj. 

where K{y,p^) = 1 + ask{y,p^) is a rapidity and transverse momentum dependent K- 
factor. We require K to be smooth for small p^, so that the smooth shape of the shower 
and MEC approaches is preserved, and to be such that 

dydp^ 

For example, we could require 

d(T< 



, Hda^LO+PS da(NLO) 



i^(y,p?)=if(y) = ^j^, (17) 



dy 



which, using the fact that 



dPT — nr = —5 — ' (1^) 

dydpT dy 

satisfies Equation [161 With this choice, the large transverse momentum tail of the distri- 
bution differs from the pure NLO result by terms suppressed by a further power of as, i.e. 
terms of NNLO order. But this does not spoil NLO accuracy, and it is therefore allowed. 
Alternatively, the K factor may also have some transverse momentum dependence, such 
that the corresponding curve in Figure [5] would be above the MEC result for px -^h, but 
would approach it for larger transverse momenta. This NNLO ambiguity in the definition 
of an NLO+PS generator is unavoidable, and it should thus be kept in mind that different 
NLO+PS generators, besides differing because of different rcnormalization and factorization 
scale choices, and different shower algorithms, may also differ because of this. 



5 MC@NLO 

In the parton shower Monte Carlo approach, one starts from the Born cross section B{^b) 
and adds higher-order corrections in the shower approximation. If one started instead from 
the NLO cross section, there would be double counting because the showers would add terms 
that are already present in the NLO result. The aim of the MC@NLO scheme is to remove 
from the NLO expressions those terms that will be generated by the parton showers. This 
is achieved by modifying the subtraction terms of the NLO calculation. The method was 
worked out in detail in [3J and applied there to vector boson pair production, then extended 
to heavy quark pair production in [TT]. In [T^ the decay angular correlations in these 
processes were added. Single top quark production processes were implemented in [131 1141 
[15] . The present version [161 HZ] also includes the single and associated production of Higgs 
and vector bosons. All these processes are implemented for hadron-hadron collisions. Heavy 
quark photoproduction has been considered in The aMC@NLO project, developing a 
fully automated NLO event generator based on the MC@NLO scheme, with loop corrections 
from MadLoop [TOj using the OPP method , has provided results on a range of more 

complex processes of interest [22l [23l [24]. A variant of the MC@NLO approach has been 
implemented within the Sherpa event generator I26j . 

Referring to Section [3] one can see that two types of NLO terms are generated by parton 
showering and need to be removed from the NLO calculation. 

1. A resolvable real emission gives rise to a positive term of the form in Equation [S] which 
has to be subtracted from the corresponding term in Equation [1] That is, one must 

replace i?($B, $rad) by i?($B, $rad) - i?(^^'($B, $rad). 
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2. In addition, expansion of the Sudakov form factor for no resolvable emission, As (Q, Qo)in 
Equation[71 to first order gives a negative term which is the integral of —R^-^'^') ((i)^ , <i)].ad) 
over the real emission phase space ^rad- 

Thus in place of Equation [T] one should start the parton shower Monte Carlo from the 
modified NLO cross section 

+ (i?($B, $rad) - i?(^^H*B, $rad))d$B d$,ad , (19) 

where the extra terms are understood to be summed over all coloured external lines of the 
Born process. By comparison with Equations [3] and 01 one can see that this amounts to a 
modification of the counterterms C($r) in the subtraction method. It is clear that, upon 
integrating over the real emission phase space $iad, Equation [19] reproduces the equiva- 
lent NLO cross section in Equation [T] However, it should be emphasised that dcTmod as 
it stands does not represent a physical differential cross section: for example, the distri- 
bution of real emissions is completely wrong and must be supplemented by the emissions 
generated by a parton shower generator. Furthermore Equation [12] is specific to a particular 
generator, represented by the label (MC). For each shower generator, one must calculate 
analytically exactly what the program does at relative order as and modify the NLO coun- 
terterms accordingly. At present there are complete MC@NLO versions [TC] [37] for the 
HERWIG UHl US] and HERWiG-h-h [301 [SI] event generators, while aMC@NLO is also 
available [221 [21] for PYTHIA [22] with virtuality-ordered parton showering. 

In the analytical calculation of the modified counterterms, to avoid the useless generation 
of unresolvable emissions and cancelling virtual corrections, the resolution cutoff Qo should 
be set to zero. This has the effect of cancelling the divergences of the pure NLO counterterms, 
provided the Monte Carlo distribution is accurate in the soft and collinear regions. 

In collinear regions, factorization of the cross section according to Equation [5] is exact and 
so the cancellation of collinear singularities is guaranteed. In soft, non-collinear regions, 
the azimuthal averaging of the splitting functions and the angular-ordering treatment of 
coherence used in the Monte Carlo mean that the cancellation of soft divergences occurs 
after integration over angles and not point- by-point. Smoothing functions are therefore 
introduced to match the MC and NLO expressions in these regions. This is not a problem 
as physical observables are insensitive to the angles of soft gluon emission. In addition, the 
factorization of colour structure in Equation [6] does not hold in general in soft, non-collinear 
regions, so there may be colour-suppressed contributions with non-cancelling sub-leading 
divergences, which give rise to power-suppressed corrections after smoothing and do not 
affect the NLO accuracy of predictions. A proof of this point is given in Appendix [XJ For 
technical details, see the Appendices of refs. [51 [TT]. 

The cancellation of divergences makes it possible to prepare separate samples of Born- 
like and real emission configurations for processing by the shower generator, known as S 
events and H events respectively, each with finite weights equal to the relevant coefficient in 
Equation [TO! This is in contrast to a pure NLO calculation, where the positive divergences 
of the real emission distribution are cancelled by negatively divergent contributions in the 
Born phase space. The S and H events can be unweighted, if desired, by accepting them 
in proportion to their weights before showering. However, these weights are not guaranteed 
to be positive, so some configurations generate counter-events that have to be subtracted 
rather than added in histograms. As can be seen from Equation [TH] provided the Monte 
Carlo is close to the true emission distribution in the important regions, the counter-events 
tend to be a small fraction of H events that serve to correct predictions to the NLO level. 

In Section[4]we claimed that a NLO+PS generator should spread the NLO K factor over 
a certain region of the transverse momentum distribution of the hardest radiation. If we 
dropped the second and third term in the first line of Equation [T^ the MC@NLO result 
would approach that of a MEC accurate Monte Carlo, with the S events being just the shower 
events, and the 11 events correcting the shower hardest emission at large transverse momenta, 
in order to match the tree- level result. The integrated cross section would be near the Born 
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cross section. In fact, the second line of Equation [19] gives only a minor contribution to 
the integral, since it is non- vanishing only in the large transverse momentum region where 
the coUinear approximation differs substantially from the tree-level matrix element. By 
including the full round bracket on the first line of Equation [THl we increase only the S 
contribution by a factor depending on the Born kinematics $b- We thus see that in the 
MC@NLO case the K factor acts uniformly over the px region that is dominated by § 
events, which typically extends out to values of the order of the hard scale of the process 
in question, while the region of harder emissions is not affected by it. 

The practical implementation of the MC@NLO method proceeds as follows. The NLO 
part first performs the integrations necessary to determine the weights of the § events and 
then generates the S and H events (i.e. Born-like and Born-|-one-parton configurations, 
respectively) to serve as the starting-points of the corresponding MC generator. These may 
remain weighted or can be unweighted as described above. For the purposes of parton 
showering and hadronization, each event has to be assigned a unique colour flow, which 
corresponds to the large- A^c limit of QCD, where Nc is the number of colours. Therefore a 
colour flow is selected according to their relative probabilities in the large- iVc limit of the 
corresponding Born or real emission matrix element. However, it should be emphasised 
that the sum of all colour flows reproduces the full NLO result, including all subleading Nc 
dependence. Each coloured external line of the event can then be processed by the shower 
generator in the normal way. In particular, there is no restriction that shower emissions 
from HI events should be softer than the extra parton emitted at the NLO level, because 
the weights of those events are computed assuming unrestricted showering. However, H 
events are not singular in the collinear and soft regions, and thus their contribution to those 
regions is phase-space suppressed, so that the cross section for producing events in which 
the shower generates radiation much harder than that generated in the H configuration is 
power suppressed. 



p'« (GeV) 




logio(p'?/GeV) ' ° ' A0W (rad) 



Figure 3: Transverse momentum (upper left plot) and rapidity (upper right) distributions 
of the top quark, and transverse momentum (lower right) and relative azimuth (lower left) 
distributions of the tt pair at the LHC (14 TeV), obtained at NLO, with HERWIG, and 
with MC@NLO. Figures from [n]. 

Figure [3] shows some MCdNLO results on top quark pair production at the LHC, 
compared with those obtained at NLO and with HERWIG. The MC@NLO and HERWIG 
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results, but not of course the NLO ones, include parton showering and hadronization. For 
ease of comparison, the HERWIG results have been rescaled by ctnlo/clO: so that ah plots 
have the same normahzation. 

For the inclusive transverse momentum and rapidity distributions of the top quark, 
higher-order corrections tend to have a global effect without causing much change of shape 
in the distributions, and consequently the results of all three approaches are similar. It 
should be noted, however, that the MC@NLO results are more stable than pure NLO with 
respect to statistical fluctuations. This is due to the absence of divergences in the weights 
within the § and IH event samples separately, rather than cancellation of oppositely divergent 
weights between the two samples. Thus it can happen that a given statistically significance 
is achieved more quickly, in spite of the extra computation involved in parton showering. 

For quantities sensitive to correlations between the top quark and antiquark, differences 
are more apparent. The MC@NLO prediction for the transverse momentum distribution 
of the pair follows the NLO result at high pT, while the HERWIG distribution, generated 
by parton showering alone, falls off rapidly at values above the top quark mass. At low 
the opposite is the case, as S events dominate and the distribution is controlled by parton 
showering, so that MC@NLO agrees with the (rescaled) Monte Carlo prediction. On the 
other hand, both NLO and HERWIG fah below the MC@NLO prediction at low values of 
the azimuthal angle between the pair, where both single hard and multiple soft emissions 
contribute. 



The basic idea in POWHEG (the acronym stands for Positive Weight Hardest Emission 
Generator) is to generate the hardest radiation first, and then feed the event to any shower 
generator for subsequent, softer radiation. In shower generators ordered in transverse mo- 
menta, the hardest emission is always the first, and in this case POWHEG simply replaces 
the hardest emission with its own, NLO accurate emission. In angular ordered showers, the 
hardest radiation may not be the first, and, as shown in ref. |3], the inclusion of so-called 
"truncated showers" , to be discussed in Section [TUl is needed to restore soft coherence in 
these cases. In POWHEG events can be produced with positive (constant) weight. Further- 
more, since the algorithm does not depend upon a particular parton shower program, the 
POWHEG output can be easily interfaced to any modern shower generator that is capable 
of handling user processes (typically those that comply with the Les Houches Interface for 
User's Processes [M]). 

The first proofs of concept of POWHEG were the implementation of Z pair production 
in hadronic collisions [22] , followed by the implementation of heavy flavour production [35] . 
The method was described in great details in ref. [5] . Drell-Yan vector boson production [37] , 
Higgs boson production via gluon fusion |38| and single-top production [39j were soon imple- 
mented. In ref. 40 , a package for the implementation of POWHEG for generic processes 
was presented, called the POWHEG BOX, that allows one, given the NLO matrix elements 
for a process, to build its POWHEG implementation automatically. Several processes have 
been implemented within this framework [4Tl|42|43l|4l|45l|46l|47l|48l|49]. Notably, in 
refs. [50", '51] m [53], the POWHEG BOX was used in conjunction with the HELAC-NLO 
package [54] for the computation of processes of considerable complexity. 

Independent POWHEG efforts have also being pursued by the Herwig+-|- [50] and 
by the Sherpa [5S] collaborations. Besides having implemented several processes in the 
POWHEG framework (SS] [57| [Si |M1 ED] , the HERWiG-h-h team has also developed the 
implementation of truncated showers, needed to recover soft coherence when interfacing 
POWHEG with an angular ordered parton shower generator. The Sherpa collaboration 
has developed a partially automated procedure for the implementation of POWHEG in |61j. 




As already discussed in Section [J] the MEG approximation in the shower formalism 
reproduces the NLO cross section at large transverse momenta. Multiplying the MEG cross 
section for the hardest emission hy a K factor that is a function of the underlying Born 
kinematics (the Higgs rapidity in our Higgs example), it is possible to achieve NLO accuracy 
for inclusive quantities. To be more precise, let us write the shower algorithm formula for 
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the generation of the largest radiation in the shower-MEC approximation 



dydpT 
A(pt) 



exp 



R 

5(Pt)A(Qo) + A(pT)^d$rad 

JD 

/R 
— d$rad6'(pT(*rad) "Pt) 



(20) 



We now claim that formula [501 achieves NLO accuracy if we replace the prefactor 

d$BS d^sB , B = B + V + j i?d<I>i.ad . (21) 

Before proving this result, we first introduce the POWHEG procedure in its full generality. 
One splits the real cross section into two components 



R — R? + R^ 1 



(22) 



where R^ is regular in the small region, and R^ embodies all the singularities. The 
POWHEG formula for the generation of the hardest radiation is 



do- = d^^B'^ 



R^ 

As(Qo) + As(pT)-^d$,.ad 



i?^d$r 



B^=B + V + J R" d$,.ad , As (pt) = exp 
A simple way to achieve the R^-R^ separation is to choose 



- / — d$,.ad6'(pT($rad) ~Pt) 



Pt 



R' 



Pt 



/l2 



Pt 



(23) 
(24) 

(25) 



With this separation, taking h oo, R^ vanishes, and one recovers the case of Equation [20l 
with the replacement Equation 1211 

In POWHEG, the hardest radiation is generated using Equation!^ This event is then 
fed to a shower Monte Carlo, which is required not to generate any radiation with larger 
transverse momentum. Because of the presence of the cutoff Qq, a small fraction of events 
without radiation will also be generated. In this case, also the shower will not generate 
radiation. 

We now demonstrate that the POWHEG formula. Equation [221 yields NLO accuracy 
when applied to IR-safe observables. We call O such an observable. First of all, we notice 
that further showering, beyond the hardest radiation, must affect the shape variable only 
to a subleading level, and it is thus enough to prove NLO accuracy for the POWHEG 
generated event. We have 



(O) 



1 
a 



0(<fB)As(Qo) 



As(pT)^0($R)d<I',.ad 



0($R)i?Pd$R 



(26) 



From Equation 1261 by adding and subtracting the same quantity, we obtain immediately 



{O) 



0($b)As(Qo)+ / As(pT)^0($B)d$,ad 



R 



d$RB=As(pT)-g- (0($r) - 0(*b)) + J 0($R)i?M$Rj . 

d$BS'50($B) + / d$Ri?^ (0(*r) - 0($b)) + / 0($R)i?^d$R 



(27) 
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By unitarity (see Equation [14]), the expression in the square bracket in the first line of Equa- 
tion[27]has been replaced by ^^(^b) in the last equality. Furthermore, another simplification 
has been applied in the last equality, which follows from the fact that the first term in the 
middle line of Equation [57] is finite (since the factor ©(^r) — 0{^b) damps the singular re- 
gion) and of NLO order. One can then drop the /B ratio and the Sudakov form factor in 
this term, since they both differ from unity by higher-order terms, which lead to corrections 
of order higher than NLO. Finally, replacing B^ with its explicit expression, we get 



(O) = -<^ / d$B B + V+lR^A^ 




='rad 



0($b)+ j d^j^R"^ (0($r) - 0($b)) 



i d$B[B + V]0{'^b) + j d<^^[R^ + i?P]0($R) 



(28) 



which concludes our proof. In the last line we assume that some sort of infrared regulator is 
applied to the virtual term and to the integral of the real term, so that infrared cancellation 
takes place as usual for infrared insensitive quantities. As we will see in the following, further 
showering is forced to be softer than the hardest radiation in POWHEG, so that IR-safe 
observables are affected by it only at a subleading level. 

It is interesting to ask ourselves what happens if we take the limit /i 0. It turns 
out that the B^ function will develop a negative infinity, since the integration of R^ has 
a positive soft-collinear singularities that cancels an analogous negative singularity in the 
virtual term. By gradually reducing i?^, the negative infinity in V will be gradually exposed. 
It is also clear that the transverse momentum distribution will approach more and more the 
NLO result, and the proof of NLO accuracy that we have just given guarantees that the 
NLO result is recovered. From this discussion, it becomes clear that h should be such that 
the nice shape of the pt distribution in the region of the resummation should not be spoiled, 
i.e. h cannot be taken too small. 

The POWHEG generated event is fed to a general-purpose shower Monte Carlo (SMC) 
such as HERWIG or PYTHIA for the generation of the rest of the parton shower and 
hadronization. Care must be taken, however, that the SMC should not generate radiation 
with transverse momenta larger than that of the initial parton. An SMC complying with the 
Les Houches Interface for User Processes (LHIUP from now on) can be easily instructed to 
do so, just by setting the common block variable scalup to the transverse momentum of the 
POWHEG generated radiation. When we talk about the transverse momentum generated 
by the SMC, we mean here either the transverse momentum with respect to the beam axis, 
for initial state radiation partons, or the transverse momentum with respect to a final state 
parton, if the emission is from a final state parton. Observe that, in the case of no radiation 
(i.e. the event consists of the Higgs boson alone), scalup has to be set to such a small value 
that further radiation from the SMC is essentially prohibited. In the following we will refer 
to the event generated by POWHEG as the Les Houches event (LHE), precisely because 
POWHEG passes the event to the SMC by using the LHIUP. 

We will not try to illustrate here how the POWHEG event is generated in practice. 
The interested reader can look at the first POWHEG implementation [3S], or at the gen- 
eral POWHEG papers HO] for details. We only note here that, besides generating the 
event kinematics, one must also assign colour to the final state partons, according to the 
LHIUP convention. The only requirement for the correctness of the POWHEG procedure 
is that colour assignment should be at least as accurate as in a shower Monte Carlo. In the 
POWHEG BOX, for example, a process-specific routine assigns the colour to the underly- 
ing Born configuration. This is generally done in the large limit, but subleading colour 
connections may also be included, if desired. The POWHEG BOX then assigns colour to 
the generated hardest emission configuration, assuming the large- A^c limit colour assignment 
in collinear radiation that is used in parton shower algorithms. 
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7 Use of parton density functions in NLO+PS genera- 
tors 



In order to claim NLO accuracy, an NLO+PS generator must make some use of NLO 
accurate parton densitites. In both MC@NLO and POWHEG the generation process 
begins with the calculation of an inclusive cross section, which is given by the two terms 
in round bracket in Equation [19] for MC@NLO, and by the calculation of B and (see 
Equation [231) in POWHEG. The B function in POWHEG, as weU as the term in round 
bracket in the first line of Equation[T2]in MC@NLO, must be computed using NLO accurate 
parton density functions (PDFs). Strictly speaking, this is not needed for the term in 
POWHEG and for the H events in MC@NLO (i.e. for the second line of Equation [TO)) , 
although in practice NLO PDFs are also alway used there. Parton densities also play a role in 
the generation of radiation via backward evolution in MC@NLO. In POWHEG they play 
a role in the generation of the hardest radiation for initial-state radiation, and subsequently 
in the SMC that takes care of the rest of the radiation. In these steps, NLO PDFs are 
not needed. The key observation needed to justify the above statements is that NLO PDFs 
are needed when they multiply a contribution to the cross section that is of leading order. 
For example, in Higgs production, NLO PDFs are needed in the contributions to the cross 
section that start at order ag. This is because NLO and LO PDFs differ by terms of order 
as, which multiplied by a term of order ag lead a term of NLO accuracy. In the NLO+PS 
generators this is only the case for the § and S cross sections in MC@NLO and POWHEG, 
respectively. In MC@NLO it is in general preferred to run the shower stage using the MC 
internal PDF's, although a user can optionally use the same NLO PDFs used in the NLO 
calculation. In POWHEG, the generation of the hardest radiation is usually performed 
using the NLO PDFs, although other choices are possible. 

8 Spin correlations in decays 

Basic processes that include narrow resonance decays can be dealt with in a standard way in 
both MCdNLO and POWHEG, considering the decay as part of the Born process. Thus, 
for example, for Z production, one can consider qq — > Z M'~M''''- In this way, the angular 
correlation of the final state decay products with the rest of the event are accounted for 
correctly. 

In the case of the production of spinless particles, as in Higgs production, it is always 
more convenient to treat the decay at a later stage, since the production formulae are 
simpler for the undecayed object, and the whole process implementation becomes simpler. 
One first generates the resonance, and then replaces it with its decay products, distributed 
isotropically in its rest frame. The same procedure can also be applied to the production 
and decay of resonances with spin, but, in doing so, one loses spin correlations. 

In ref. [12] , a method for the inclusion of spin correlations in the decay of resonances 
in Monte Carlos, applicable also to NLO+PS implementations, has been introduced. In 
essence, the method works as follows. The MG@NLO framework generates an event to be 
passed to the shower Monte Carlo, which is simply the partonic Born configuration in the 
case of § events, or the Born configuration plus the radiation of one parton in H events. 
The idea is essentially to complete the partonic event by letting the heavy resonance decay 
and weighting the angular distribution of its decay products according to the corresponding 
tree-level matrix element. This procedure certainly captures correctly the Born level angular 
distribution, and, for hard radiation, also the correct correlation when an associated hard 
jet is produced. The same procedure can be applied to the POWHEG case. In POWHEG, 
most events include the hardest radiation. In this case the correlation is retained even if the 
radiated parton is near the coUinear region, or very soft. It should be kept in mind, however, 
that the procedure of ref. p[2 does not include the full NLO correction to correlations, since 
the spin dependence of the virtual corrections is not fully accounted for. 

In final state resonance decays into coloured particles, like qq — ^ Z — > qq, or top pro- 
duction and decays, the same method can be, and in fact is, used. Chain decays, such as 
t — !> Wb qq, are also handled. However, next-to-leading corrections to the decay process 
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are not included. Thus, for example, in the implementation of heavy flavour production in 
MC@NLO and POWHEG, the top decay is only handled at the leading-order level, and 
only at this level are spin correlations accounted for correctly. 



9 MC@NLO and POWHEG for complex processes 



The discussion given so far has focussed upon a relatively simple process, where only initial- 
state coUinear radiation has to be dealt with. In general, final-state massless partons may 
also be present, and should be properly considered. Typically, processes like single top pro- 
duction, or Higgs production via vector boson fusion, have final-state light partons that can 
emit coUinear radiation. In processes of associated jet production, like vector boson plus 
one jet production, or dijet production, one has the further complication that the Born cross 
section itself is divergent, unless some transverse momentum cut on the jet kinematics is 
imposed. All these issue add practical complexity to the NLO-I-PS implementation. How- 
ever, no conceptual complexity arises: the same method used in the simplest cases can be 
applied in the most complex ones. The key observation is that the real cross section can 
be split into a sum of contributions with specific flavour structure that are singular in one 
coUinear region only. One writes 



where each a labels a particular flavour structure and a singular region of the real amplitude, 
and the i?" is required to be singular only in the corresponding region. Considering, for 
example, Z-fjet production, the real-emission flavour structure qq — )■ Zgg has the following 
singular configurations: the first emitted gluon may have vanishing transverse momentum 
(i.e. it is soft or coUinear to either of the initial-state partons), the second emitted gluon may 
have vanishing transverse momentum, or the two emitted gluons may have small relative 
angle. The full cross section for the qq Zgg process may be written as 



with q > 2, where is the transverse momentum of the first gluon, is the transverse 
momentum of the second gluon, and Pj, is the relative transverse momentum of the two 
final state gluons. It is clear that i?i is singular only when the first gluon is coUinear to 
either initial-state parton, i?2 is singular only when the second gluon is coUinear to either 
initial-state parton, and R3 is singular only if the two final-state gluons are coUinear to each 
other. Notice also that the soft singularities are partitioned among the different coUinear 
regions, but that in all cases nothing is omitted, since the sum of the three contribution 
yields the total real amplitude by construction. Both POWHEG and MC@NLO adopt a 
similar separation of the coUinear regions. In the computation of the § events in MC@NLO, 
the third term on the first line of Equation \W\ should include all contributions arising from 
the shower, i.e. initial-state radiation and final-state radiation. The term on the second 
line of Equation [19] is instead split into all singular components of R, and each component 
is accompanied by the corresponding MC component. Thus, for example, the final-state 
radiation contribution will have the structure 



Carlo implementation of the final-state splitting process. Notice that the NLO and MC 
terms in Eg nation 13 1 1 must be expressed in terms of the same variables, including appropriate 
Jacobian factors. 

In POWHEG, in the case of multiple singular regions, one must first of all define the 
appropriate B function as 




(29) 



a 




(30) 





(32) 
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Here fb represents a specific Born flavour configuration, and the notation {a|/fc} stands for 
all a such that i?^ has ft as flavour of the underlying Born configuration. Notice that the 
mapping of $b and $rad into the full real phase space also depends upon the particular 
singular contribution that we are considering, and the superscript a on $rad is there to 
remind us of this fact. We can think of the -B§ as the inclusive NLO cross section at fixed 

Jb 

underlying Born flavour ft, and kinematics $b- Of course, the definition of the underlying 
Born structure also depends upon the separation of the real cross section R into the singular 
components R^, and upon the mapping that is used in a particular singular region for 
constructing the real phase space out of the underlying Born phase space and the radiation 
variables. A detailed description of the decomposition Equation [^S] and of the mappings is 
given in ref. [S]. 

The generation of radiation in POWEEEG, once the underlying Born kinematics and 
flavour structure have been generated, is obtained using as Sudakov form factor the expres- 
sion 

Again, the sum in the exponent is over all a such that R^ has as the flavour of the 
underlying Born configuration, and the integral is performed at fixed underlying Born con- 
figuration. The function fc^"'' is the radiation transverse momentum, i.e. the transverse 
momentum of the radiated parton in initial-state radiation, and the transverse momentum 
of the radiated parton relative to a final-state parton in final-state radiation, and for this 
reason it depends upon the singular region under consideration. The Sudakov form factor in 
Equation 133] is the probability that no radiation, either initial- or final-state, from any of the 
final-state partons, has fcx larger than px- In practice, in order to generate the radiation, the 
full Sudakov form factor is decomposed in the product of Sudakov form factors for each i?^ 
contribution. According to standard Monte Carlo methods, the generation of px according 
to Equation (331) can be performed by generating a px using each Sudakov form factor, and 
picking the largest. 

The generation of the contribution of the non-singular R^ terms does not present any 
particular problem. 



10 Truncated showers 



In angular-ordered showers, like those in HERWIG and HerwigH — h, the hardest emission 
during the shower development may not necessarily be the first in the shower chain. As seen 
in Section [31 angular ordering allows softer, larger-angle emission to take place earlier in 
the shower development, in order to account for coherent soft gluon emission from a bunch 
of coUinear partons. These soft, large-angle emissions see the bunch of coUinear emitting 
partons as a single coherent colour source, its colour being that of the parent parton of the 
bunch. This is why the shower can treat these emissions as if they were coming from a 
single parton. In ref. i4|, it was shown that even in the case of angular-ordered showers, the 
hardest parton emission is described, up to subleading corrections, by the equation 



dcr 



A(pt) 



Bd$B 



exp 



A(Qo) + A(pT)^^d$,ad 



B 



B 



(34) 



which has the same form as Equation [T31 It was thus suggested that the Herwig shower was 
equivalent to a shower initiated by Equation I34p . which we call in the following the initial 
process, with the remaining radiation provided as follows: 

I A shower initiated at an angle determined by the underlying Born configuration of the 
initial process, stopping at an angle equal to that of the emission in the initial process 
, vetoing radiation harder (i.e. with larger transverse momentum) than the initial- 
process radiation. This shower was called a vetoed truncated shower in ref. 
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II All partons produced in the initial process are allowed to radiate with standard angular- 
order initial conditions, but radiation with transverse momentum larger than that of 
the initial process is vetoed. 

The vetoed truncated shower described in item I is needed in order to supply the soft, 
coherent radiation that can be emitted at large angles by the splitting partons associated 
with the hardest emission. An angular-ordered shower algorithm generates this radiation 
first, so the hardest emission may take place later in the shower. If instead the shower is 
started with the hardest emission, this coherent radiation must be added explicitly. 

In ref. [3] it was proposed to implement a positive weights NLO-I-PS generator by gener- 
ating the initial process with the POWHEG method, and to complete the shower according 
to the prescriptions I and II. It was also pointed out that vetoed-truncated showers are 
not really specific requirements of POWHEG, but that they also appear naturally in the 
CKKW matrix-element to shower matching (ME-I-PS). In ref. [62 , an implementation of 
CKKW matching using truncated showers has been proposed in the HerwigH — h framework. 
An implementation of truncated showers in the Sherpa framework has been proposed in 
[55] . However, since Sherpa uses a transverse momentum ordered dipole shower, one can 
no longer claim that truncated showers are needed there to restore the coherence of soft 
emission. 

11 POWHEG and MC@NLO comparisons 

In the present section we compare and contrast the MC@NLO and POWHEG methods. As 
discussed in Section [TUJ it was shown in ref. [3] that the Herwig shower could be replaced by 
a shower where the hardest emission is generated first, and that the hardest emission, up to 
subleading corrections, is described by Eauation l34l It then follows that, up to their nominal 
accuracy, we can identify the MC@NLO prescription with a POWHEG implementation 
having = R^^^\ where the § and the H events in MC@NLO are identified with the S 
and F events respectively in POWHEG. This reasoning provides an alternative to the proof 
in ref. [3] of the correctness of the MC@NLO procedure by showing that it is equivalent to 
a POWHEG approach, so that the proof of NLO accuracy for inclusive observables used for 
POWHEG also applies. Furthermore, it allows us to better understand the similarities and 
differences between the two approaches. As will be made clear, they can mostly be traced 
back to differences between and /J^'^^' . 

Detailed comparisons of MCQNLO and POWHEG were carried out in refs. [551 1551 1551 
1371 139] . Remarkable agreement was found for inclusive quantities, i.e. quantities that are 
inclusive in the hardest jet being radiated. Discrepancies were found in the kinematics of the 
hardest jet, especially for processes with a very large NLO K factor, like Higgs production. 
We will focus now upon these distributions, bearing in mind that all the discrepancies reflect 
differences in the estimation of corrections that are beyond the formal accuracy of the two 
approaches. 

We begin by showing the transverse momentum distribution of the Higgs in the left plot 
of Figure 21 The POWHEG and MC@NLO predictions agree for moderate transverse 
momenta. For large transverse momenta POWHEG gives a larger results, with the upper 
POWHEG prediction (corresponding to ft. = cxd) in the last bin in the plot being above 
MC@NLO by roughly a factor of three. This difference can be understood as due to 
the part of the transverse momentum spectrum that is amplified by the NLO K factor in 
the two programs. For transverse momenta below about 100 GeV, the MC@NLO result is 
dominated by § events, and is thus equal to the pure HERWIG result multiplied by an NLO 
K-iactOT. At larger transverse momenta, the H events start to contribute, and eventually 
dominate the spectrum. The difference with respect to POWHEG is due to the fact that 
in the latter, for h = oo, one chooses R^ = R and R^ = 0. Therefore, the NLO K factor 
multiplies uniformly the whole transverse momentum distribution, amplifying it by a factor 
of order 2, even in the region of high transverse momenta, where MCdNLO is dominated 
by H events and no X-factor is present. Notice also that the POWHEG prediction with 
h = 100 GeV is below the h = oo prediction by almost a factor of 2, since with this choice 
that region is dominated by F events. A remaining factor of about 1.6 arises from the fact 
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Figure 4: Transverse momentum distribution of the Higgs (left plot) and distribution of 
the rapidity distance between the Higgs and the hardest jet, obtained with HERWIG, with 
MC@NLO and with MCONLO with the the NLO corrections to the S events switched 
off. For comparison, the POWHEG+PYTHIA resuh is also shown for the h ~ oo and 
h = 100 GeV value of the parameter of eq. [511 



that MC@NLO uses as default cent ral value f or the factorization and renormalization scales 
the Higgs transverse mass mJj! = + , while in POWHEG the Higgs mass is used. 

It is easily seen that, for the given choice of center-of-mass energy and Higgs mass, the gluon 
density is quite insensitive to the factorization scale, and the whole remaining effect is due to 
the different choices of renormalization scale in the Q;g(/xp{ ) dependence of the cross section 
at large transverse momentum. 

We report in the figure also MC@NLO with no NLO terms included in the S events. 
In this way, one gets essentially the MEC corrected spectrum, with the hard tail matching 
the fixed-order NLO result, but the total rate nearly equal to the LO result. We can see 
that this result is closer in shape to the POWHEG result, the difference being due to the 
different renormalization scale choice. 

The right plot of Equation H] shows the distribution in the rapidity difference between the 
hardest jet and the Higgs boson for a jet pT cut of 20 GeV. In this case, the pure HERWIG 
result has a dip at zero rapidity difference, a feature due to the shower approximation, since 
radiation is built from the independent radiations of the two incoming partons. This feature 
is somewhat sensitive to the shower generator (see e.g. [HI])- From the figure it is apparent 
that in the MC@NLO result with no NLO correction to the § events, the inclusion of the 
H events corrects almost completely for this effect. This is not unexpected, since the HI 
events are generated with a distribution equal to the difference between the real one-parton 
emission and the HERWIG approximation to it. In the full MC@NLO result, however, § 
events are amplified by the large K factor, which amplifies the dip in the pt region where 
S events contribute significantly, and the H events, which are not multiplied by the same K 
factor, can no longer compensate for it. Notice however that this effect is in fact of order ag, 
i.e. it is beyond the NLO approximation. It was first noticed in ref. |64j in the framework 
of tt production. An explanation of the effect along the lines reported here has been given 
in refs. [3S1 [51]. Even if somewhat extreme, the effect is within the uncertainties due to 
unknown NNLO terms. It unfortunately affects the shape in an unphysical way, since no 
behaviour of this kind is observed in the 0{ag) calculation of the Higgs plus one jet cross 
section. It reminds us, however, that also the shape of distributions is subject to the same 
NNLO uncertainties as those affecting the total cross section. 



12 Uncertainties in NLO+PS 

Several uncertainties affect NLO+PS generators. Some of them are specific to the shower 
and hadronization parts, and others are related to the parameterization of the parton den- 
sities and to the chosen value for the strong coupling constant. Here we want to focus upon 
the uncertainties that are characteristic of the QCD NLO calculation. We know that fac- 
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torization and renormalization scale uncertainties will play an important role in estimating 
the effect of higher order terms. There are, however, other sources of uncertainty that are 
specific to the NLO+PS algorithms, and in fact to all NLO calculations that are improved 
with the leading log resummation of soft gluons effects. As discussed in Section HI the dis- 
tribution of the transverse momentum of the radiation has a singularity at zero momentum, 
and the value of the NLO corrections is determined by the interplay of the soft divergences 
in the virtual corrections and of the singularity of the real cross section at zero transverse 
momenta. On the other hand, NLO+PS generators, as well as resummation formulae, turn 
this singularity into a smooth curve. NLO effects are thus spread over a wider region, and 
a further scale will implicitly appear to delimit this range. We have seen in the previous 
section that this implicit scale is the cutoff scale for the MC radiation in MC@NLO (the 
Higgs mass in the Higgs production example), while it can be up to the kinematic limit for 
the transverse momentum of radiation in POWHEG. 

NLO uncertainties for quantities that are inclusive in the hardest radiation can safely be 
estimated using standard scale variation. On the other hand, for quantities that are sensitive 
to the hardest radiation some care is needed. In fact, for S or S events in MC@NLO or 
POWHEG respectively, the shape of the transverse momentum distribution of the hardest 
radiation is unaffected by scale variation. In practice, scale variation in the square bracket 
on the r.h.s. of Equation!^ is never performed, since by doing it one easily spoils the NLL 
accuracy of the Sudakov form factor. Similarly, a scale variation of the corresponding term 
in MC@NLO is never performed, since it can only be achieved by changing the scale in the 
Monte Carlo event generator that is being used. Thus, the scale variation in the S events in 
POWHEG, and in the § events in MC@NLO, only affects the B prefactor. This implies 
that the effect of scale variation on S (or §) events is of relative order ag, while that of F 
(or H) events is of relative order as. This point can be better illustrated as follows. The B 
prefactor is of order q;| at the Born level, and it includes NLO corrections, of order ag. Its 
scale variation must therefore be of order ag, i.e. beyond its nominal accuracy. Therefore 
the relative scale variation 5B/B is of order ag. On the other hand, the F(H) term is of 
order ag, and its scale variation is of order ag, so its relative scale variation is of order as. 
Thus, the larger the contribution to the transverse momentum distribution coming from S 
(or §) events, the smaller its relative scale dependence will be. In the following, we will focus 
on this problem, taking as an example the transverse momentum distribution of the Higgs 
boson. This can also be computed (using the HqT program [65l [66l EZl [68] ) at a matched 
NNLL+NNLO accuracy, and can thus serve as a benchmark for assessing how realistic are 
the NNLO terms introduced by the NLO+PS approaches. We will include in our calculation 
the perturbative uncertainty determined by scale variation. The studies reported here are 
taken from ref. 691. 

The plots that follow have all been obtained with the following settings. We have used the 
MSTW2008 NNLO central PDF set [70] for all curves. This is because HqT requires NNLO 
parton densities, and because we want to focus upon differences that have to do with the 
calculation itself, rather than the PDFs. The HqT result has been obtained by running the 
program with full NNLL+NNLO accuracy, using the "switched" result. The resummation 
scale Q in HqT has been set to AfH/2. The bands shown in the figures for HqT, POWHEG 
and MC@NLO are all obtained by varying the factorization and renormalization scales 
by a factor of two above and below their default central value, imposing the constraint 
0.5 < /ip/z^R < 2. In the MC@NLO case, the independent scale variation is not performed, 
and only equal scales are used. 

It is instructive to analyze the difference between MC@NLO and POWHEG at their 
default value of parameters. This is illustrated in Figure [5] As shown earlier, a large 
difference in the high transverse momentum tail is visible. As can be seen immediately, 
when the uncertainty bands due to scale variations are included, the differences between 
the two NLO+PS methods, and their differences with respect to the HqT result, become 
less pronounced. The observation made earlier on the insensitivity of the shape of the pt 
distribution contributed by S (or S) events is also visible, with the POWHEG shape being 
constant for all transverse momenta, and the MC@NLO shape being roughly constant in 
the region dominated by H events. It is also seen that the scale uncertainty in POWHEG 
corresponds to a uniform factor in the whole range of transverse momentum being considered. 
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Figure 5: The transverse momentum spectrum of the Higgs in MC@NLO (upper) and in 
POWHEG+PYTHIA (lower) compared to the HqT result. In the lower insert, the same 
results normalized to the HqT central value are shown. 



In order to give a more realistic assessment of the uncertainties in POWHEG, one can 
also exploit the freedom in the separation R — + R^ . In the case of Higgs production, it 
is found that by choosing h — Mii/1.2, the POWHEG result closely matches in shape that 
of HqT. It is similarly found that MC@NLO better matches the HqT output if Mh, rather 
than Tot, is used as central scale for scale variations. These results are shown in FigurelHl We 
see that now the large differences between MC@NLO and POWHEG are mostly removed, 
since both adopt a central scale equal to Mh, and both adopt a similar separation of S (or 
§) and F (or H) events. Both generators, furthermore, display a reasonable scale variation 
in the high regime, while the scale variation at moderate values of transverse momenta 
(around 100 GeV) seems to be comparable to that of HqT, rather than being larger. 

In MCdNLO, we have seen that the shape of the transverse momentum spectrum 
at moderate pt exhibits only a mild dependence upon scale variation. This shape is in 
fact determined by S events, and thus depends only upon the shower Monte Carlo that is 
being used, which is HERWIG in the present case. We may expect significant changes 
in shape if other Monte Carlos are used. In Figure [7] we display the Higgs pr spectrum 
using MCdNLO with the virtuality-ordered version of PYTHIA [32] ■ We do indeed see 
a considerable difference in the spectrum at small transverse momenta. The discrepancy 
with HqT at small transverse momenta is purely due to the fact that the virtuality-ordered 
version of PYTHIA does not match well with HqT at small transverse momenta. 



13 NLO+PS versus ME+PS matching 

Matching tree-level matrix elements and parton shower generators (ME-I-PS) allows for the 
generation of samples where a basic process is accompanied by a fairly large number of 
associated jets. The ME-t-PS method was first formulated in ref. [I] (CKKW), and several 
variants have appeared since (for a summary of the various implementations see refs. pil71]V 
As a representative example, taking W boson production as the basic process, the method 
allows one to construct a sample of W with an arbitrary number of associated jets, where 
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Figure 6: Uncertainty bands for the transverse momentum spectrum of the Higgs bo- 
son at LHC, 7 TeV, for a higgs mass ttih = 120 GeV. On the upper pfots, the 
MC@NLO+HERWIG result obtained using the non-dcfauh value of the reference scale 
equal to A/h- On the lower plots, the POWHEG+PYTHIA output, using the non-default 
+ separation. The uncertainty bands are obtained by changing /xr and /ip by a 
factor of two above and below the central value, taken equal to Mjj, with the restriction 
0.5 < ^r/mf < 2. 



the distributions of the first n jets {n being limited by the computing power available) are 
computed with tree-level accuracy, and the remaining ones are generated in the collinear 
approximation. The ME-fPS method does not achieve NLO accuracy for any quantity. 
However, it cannot be claimed that the NLO-I-PS methods are the NLO extension of the 
ME+PS method. In the example of W production, an NLO-I-PS implementation of W 
production can be used to produce a sample such that inclusive W distributions, such as 
the W rapidity distribution, are accurate at the NLO level. On the other hand, -t- 1 jet 
production is described by this sample at order as , i.e. with the same level of accuracy as 
an ME-I-PS implementation. Furthermore, W + 2 oi more jets is generated in the NLO-I-PS 
only with shower accuracy, i.e. in the collinear approximation, while in the ME-I-PS sample 
it is generated with tree-level accuracy. 

In ref. [64], the ALPGEN ME-I-PS generator was presented and a thorough comparison 
with MC@NLO was performed for top production. Figure [8] shows the main features of the 
comparison. Good agreement is found for inclusive quantities like the top quark transverse 
momentum distribution, provided that ALPGEN (which does not include NLO corrections) 
is rescaled hy a. K factor. Once the K factor is accounted for, the zero- and one-jet cross 
sections agree reasonably in the two approaches. At larger jet multiplicities we expect and 
see discrepancies between the two approaches, ALPGEN being more reliable in this case as 
it does not use the collinear approximation for the production of extra jets. 



14 Outlook and further developments 

The present NLO-I-PS implementations leave much room for improvement, and several pa- 
pers in the literature propose new approaches. First of all, it would be desirable to merge the 
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Figure 7: Uncertainty bands for the transverse momentum spectrum of the Higgs bo- 
son at LHC, 7 TeV, for a higgs mass rriH = 120 GeV. On the upper pfots, the 
MC@NLO+PYTHIA resuh obtained using the non-defauh vahie of the reference scale 
equal to Afn- The bare PYTHIA resuh rescaled by a iiT- factor is also shown. 
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Figure 8: Comparison of MC@NLO and ALPGEN in top production, for a 14 TeV LHC, 
from ref. [M]. On the left, the transverse momentum of the top quark. On the right, the 
jet multiplicity. 



NLO+PS and ME+PS methods, in such a way that higher jet multiplicities are described at 
tree-level accuracy while inclusive observables maintain NLO accuracy. A further goal is the 
full extension of the ME-fPS method to NLO, and several proposals in this direction have 
appeared in the literature [751 [731 [7H [751 [73] . Fixed-order NNLO calculations have become 
available for some collider processes, and their implementation in a shower framework would 
be welcome. Finally, a full extension of the shower algorithm to NLO, i.e. including NLO 
splitting kernels is being pursued [771 [ZH] • 

Besides pursuing new approaches, one can also investigate to what extent some of these 
objectives can be approached by simply merging event samples obtained with available tools. 
In ref. [79], a recipe for merging a POWHEG together with a MadGraph ME+PS sample 
is given for the cases of W and tt production, and in ref. |80j a practical recipe is presented 
for merging the Z and Z + 1-jet POWHEG samples. 
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Appendix 



A Smoothing procedure in MC@NLO 



In this appendix we demonstrate that the smoothing procedure used in MC@NLO to cure an 
imperfect cancellation of soft divergences between R and i?'^*-^) has only power-suppressed 
effects on infrared-safe observables. Although this is a technical point, we include it here 
just as an example of how the method used to prove NLO accuracy in Section |5] can also be 
applied to other important issues. First of all, in Section[TO]we have shown that MC@NLO 
is equivalent to a POWHEG generator with = R'-^'^\ and thus the proof of NLO 
accuracy given in Section [5] can also be applied for MC@NLO. 

Assume that the shower approximation for the real cross section is matched to the exact 
real cross section at a scale Qm, i-e. that for transverse momentum of the radiation below 
Qm, R^'~^ is smoothly matched to R. We write i?^^ fo^' the matched cross section. The 
proof of NLO accuracy in Section [6] goes more or less as before, and by the time we reach 
Equation it has the following form 



a{0) 



B + V+ I i?,^^d$,ad 



0($b) + J d$Ri?^C(0($fl,)-0($B)) 



d$R [R - R^,^] 0{<fR) 



(35) 



T-MC 



where we have assumed that i?^^ replaces R^'^ in the § events cross section B"'", and in 
R^ , which is the cross section for H events. Each term in Equation [35] is finite, since i?^^ 
matches R in the singular region. We can easily manipulate Equation 1351 as follows: 



a{0) 



I 



B + V 



J i?d$rad 0($b)+ J d$Ri?^C(0($ii)-0($B)) 



+ 



(36) 



where, with respect to Equation I35[ the term proportional to 0($b) subtracted from the 
last term has been added to the first, and finally 



a{0) 



B + V 



i?d$,. 



ad 



0($b) + y"d$Ri?(0($fl) -0($b)) 



y'd<i>R[i?^c-i?,'::c] [o($h)-o($b)] 



(37) 



The first line in Equation [571 corresponds to the exact NLO result, and the last line corre- 
sponds to the correction due to the smoothing procedure. Now we suppose that [i?^^ — R^'"] 
is singular in the singular region. However, [0($fl) — 0(<i>B)] kills the singularity for 
infrared-safe observables, so, the integrand in Equation [37] is finite. Furthermore it van- 
ishes if the transverse momentum of radiation is above Qm, because R^'^ = R in that 
case. Thus we have a power of Qm suppression, where the exact power depends upon how 
smoothly 0($_r) approaches 0($) near the singular region. 
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